{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "202568 005.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/202568/005.swc\n",
      "202568 011.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/202568/011.swc\n",
      "202568 013.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/202568/013.swc\n",
      "202568 017.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/202568/017.swc\n",
      "202568 021.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/202568/021.swc\n",
      "202568 022.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/202568/022.swc\n",
      "202568 009.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/202568/009.swc\n",
      "210728 001.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/001.swc\n",
      "210728 002.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/002.swc\n",
      "210728 003.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/003.swc\n",
      "210728 004.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/004.swc\n",
      "210728 005.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/005.swc\n",
      "210728 010.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/010.swc\n",
      "210728 012.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/012.swc\n",
      "210728 014.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/014.swc\n",
      "210728 019.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/019.swc\n",
      "210728 033.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/033.swc\n",
      "210728 035.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/035.swc\n",
      "210728 041.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/041.swc\n",
      "210728 042.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/042.swc\n",
      "210728 048.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210728/048.swc\n",
      "210729 001.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/001.swc\n",
      "210729 003.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/003.swc\n",
      "210729 004.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/004.swc\n",
      "210729 005.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/005.swc\n",
      "210729 006.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/006.swc\n",
      "210729 054.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/054.swc\n",
      "210729 055.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/055.swc\n",
      "210729 057.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/057.swc\n",
      "210729 059.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/059.swc\n",
      "210729 061.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/061.swc\n",
      "210729 062.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/062.swc\n",
      "210729 063.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/063.swc\n",
      "210729 064.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/064.swc\n",
      "210729 066.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/066.swc\n",
      "210729 093.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210729/093.swc\n",
      "210731 002.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210731/002.swc\n",
      "210731 007.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210731/007.swc\n",
      "210731 010.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210731/010.swc\n",
      "210731 011.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210731/011.swc\n",
      "210731 049.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/210731/049.swc\n",
      "211185 011.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211185/011.swc\n",
      "211185 014.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211185/014.swc\n",
      "211185 024.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211185/024.swc\n",
      "211185 025.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211185/025.swc\n",
      "211185 028.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211185/028.swc\n",
      "211185 031.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211185/031.swc\n",
      "211186 002.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211186/002.swc\n",
      "211187 001.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211187/001.swc\n",
      "211187 006.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211187/006.swc\n",
      "211187 020.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211187/020.swc\n",
      "211187 022.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211187/022.swc\n",
      "211188 002.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211188/002.swc\n",
      "211189 003.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211189/003.swc\n",
      "211189 005.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211189/005.swc\n",
      "211190 002.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/002.swc\n",
      "211190 003.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/003.swc\n",
      "211190 004.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/004.swc\n",
      "211190 006.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/006.swc\n",
      "211190 007.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/007.swc\n",
      "211190 008.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/008.swc\n",
      "211190 009.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/009.swc\n",
      "211190 010.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/010.swc\n",
      "211190 012.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/012.swc\n",
      "211190 013.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/013.swc\n",
      "211190 014.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/014.swc\n",
      "211190 016.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/016.swc\n",
      "211190 017.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/017.swc\n",
      "211190 019.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/019.swc\n",
      "211190 020.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/020.swc\n",
      "211190 022.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/022.swc\n",
      "211190 025.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/025.swc\n",
      "211190 027.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/027.swc\n",
      "211190 029.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/029.swc\n",
      "211190 030.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/030.swc\n",
      "211190 032.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/032.swc\n",
      "211190 034.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/034.swc\n",
      "211190 036.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/036.swc\n",
      "211190 037.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/037.swc\n",
      "211190 062.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/062.swc\n",
      "211190 067.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211190/067.swc\n",
      "211191 005.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211191/005.swc\n",
      "211192 001.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211192/001.swc\n",
      "211193 005.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211193/005.swc\n",
      "211193 006.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211193/006.swc\n",
      "211193 008.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211193/008.swc\n",
      "211193 009.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211193/009.swc\n",
      "211193 010.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211193/010.swc\n",
      "211193 012.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211193/012.swc\n",
      "211193 015.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211193/015.swc\n",
      "211466 033.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211466/033.swc\n",
      "211467 014.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211467/014.swc\n",
      "211470 020.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211470/020.swc\n",
      "211471 005.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211471/005.swc\n",
      "211985 002.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/002.swc\n",
      "211985 003.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/003.swc\n",
      "211985 004.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/004.swc\n",
      "211985 005.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/005.swc\n",
      "211985 006.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/006.swc\n",
      "211985 014.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/014.swc\n",
      "211985 035.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/035.swc\n",
      "211985 036.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/036.swc\n",
      "211985 037.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/037.swc\n",
      "211985 043.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211985/043.swc\n",
      "211988 006.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/006.swc\n",
      "211988 014.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/014.swc\n",
      "211988 016.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/016.swc\n",
      "211988 022.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/022.swc\n",
      "211988 023.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/023.swc\n",
      "211988 024.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/024.swc\n",
      "211988 030.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/030.swc\n",
      "211988 043.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/043.swc\n",
      "211988 053.swc\n",
      "exist  d:\\project\\python\\neuron-vis\\figures/../neuronVis/../resource/swc/211988/053.swc\n",
      "[4 4 4 ... 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "import sys,copy,os,inspect\n",
    "if hasattr(sys.modules[__name__], '__file__'):\n",
    "    _file_name = __file__\n",
    "else:\n",
    "    _file_name = inspect.getfile(inspect.currentframe())\n",
    "CURRENT_FILE_PATH = os.path.dirname(_file_name)\n",
    "sys.path.append(os.getcwd()+\"/../neuronVis\")\n",
    "\n",
    "import IONData,Scene,SwcLoader\n",
    "import NeuronProcess\n",
    "from sklearn.cluster import KMeans\n",
    "import numpy as np\n",
    "import joblib\n",
    "\n",
    "def GetFirstOrderEdges(edgesForCluster,parentedge):\n",
    "    for child in parentedge.children:\n",
    "        if child.order!=0:\n",
    "            edge = SwcLoader.Edge()\n",
    "            edgesForCluster.append(edge)\n",
    "            for p in child.data:\n",
    "                edge.addPoint(SwcLoader.Point([p.index,p.type,p.x,p.y,p.z,p.ratio,p.parentIndex]))\n",
    "            child.maxDepth=1\n",
    "            GetFirstOrderEdges(edgesForCluster,child)\n",
    "#In[0]\n",
    "\n",
    "\n",
    "import Visual as nv\n",
    "import json\n",
    "# neuronvis = nv.neuronVis(renderModel=1)\n",
    "\n",
    "# neuronvis.render.setBackgroundColor((0.0,0.20,0.5,1.0))\n",
    "# neuronvis.addRegion('CA1',[0.5,1.0,0.5])\n",
    "iondata = IONData.IONData()\n",
    "# f=open('../resource/pfcsubtype60.json', encoding='gbk')\n",
    "# neurons=[]\n",
    "# neurons = json.load(f)\n",
    "# neurons=iondata.getNeuronListBySampleID('211190')\n",
    "neurons= Scene.scene2List('../resource/scene/cluster5_spcd20220628_spcd.nv')\n",
    "# neuronvis.render.setLookAt((5000,0,15000),(5000,0,0),(0,-1,0))\n",
    "\n",
    "edgesForCluster=[]\n",
    "\n",
    "for neuron in neurons[0:-1]:\n",
    "    print(neuron['sampleid'], neuron['name'])\n",
    "    neuronT = iondata.getNeuronTreeByID(neuron['sampleid'], neuron['name'])\n",
    "    neuronT.dendriteHide=True\n",
    "    NeuronProcess.CalculateBranchMaxDepth(neuronT)\n",
    "    mainbranch=[]\n",
    "    if neuronT.rootAxonEdge is None:\n",
    "        continue\n",
    "    for edge in neuronT.edges:\n",
    "\n",
    "        if edge.maxDepth>neuronT.rootAxonEdge.maxDepth*0.5:\n",
    "            edge.order=0\n",
    "            mainbranch.append(edge)\n",
    "    for edge in mainbranch:\n",
    "        for child in edge.children:\n",
    "            if child.order!=0:\n",
    "                NeuronProcess.OrderChildren(neuronT,child,1)\n",
    "        pass\n",
    "    # for child in neuronT.rootAxonEdge.children:\n",
    "    # \t# if child.order!=0:\n",
    "    # \t\tNeuronProcess.OrderChildren(neuronT,child,1)\n",
    "\n",
    "    for edge in mainbranch:\n",
    "        for child in edge.children:\n",
    "            if child.order!=0:\n",
    "                GetFirstOrderEdges(edgesForCluster,child)\n",
    "\n",
    "newneuron=SwcLoader.NeuronTree()\n",
    "if len(mainbranch):\n",
    "    newneuron.edges= edgesForCluster\n",
    "    newneuron.rootAxonEdge = neuronT.rootAxonEdge\n",
    "    newneuron.root=neuronT.root\n",
    "# neuronvis.addNeuronTree(newneuron,'neuronname',depthIntensity=False)\n",
    "# neuronvis.clear()\n",
    "data = []\n",
    "for edge in edgesForCluster:\n",
    "    xyz=[]\n",
    "    edge.children=[]\n",
    "    for point in edge.data:\n",
    "        xyz.append(point.xyz)\n",
    "    # print(np.array(xyz).mean(axis=0))\n",
    "    data.append(np.array(xyz).mean(axis=0))\n",
    "cluster=5\n",
    "# print(data)\n",
    "estimator=KMeans(n_clusters=cluster,max_iter=1500)\n",
    "res=estimator.fit_predict(data)\n",
    "# 预测类别标签结果\n",
    "lable_pred=estimator.labels_\n",
    "# 各个类别的聚类中心值\n",
    "centroids=estimator.cluster_centers_\n",
    "# 聚类中心均值向量的总和\n",
    "inertia=estimator.inertia_\n",
    "print(lable_pred)\n",
    "for i in range(len(edgesForCluster)):\n",
    "    edge=edgesForCluster[i]\n",
    "    edge.maxDepth=lable_pred[i]\n",
    "neuronT.rootAxonEdge.maxDepth=4\n",
    "newneuron=SwcLoader.NeuronTree()\n",
    "# if len(edgesForCluster):\n",
    "#     newneuron.edges= edgesForCluster\n",
    "#     newneuron.rootAxonEdge = neuronT.rootAxonEdge\n",
    "#     newneuron.root=neuronT.root\n",
    "#     neuronvis.addNeuronTree(newneuron,'t',depthIntensity=True)\n",
    "# neuronvis.render.savepng('../resource/test42.png')\n",
    "# neuronvis.render.run()\n",
    "\n",
    "\n",
    "#In[1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "res,path = iondata.getFileFromServer(\"IB_MP_P_MY.nrrd\")\n",
    "import nrrd\n",
    "IB_MP_P_MYmap,header = nrrd.read(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusterLength={}\n",
    "for c in range(5):\n",
    "    clusterLength[c]={0:0,354:0,771:0,313:0,1129:0}\n",
    "for i in range(len(edgesForCluster)):\n",
    "    edge=edgesForCluster[i]\n",
    "    lenMap ={0:0,354:0,771:0,313:0,1129:0}\n",
    "    for j in range(len(edge.data)-1):\n",
    "        p = edge.data[j]\n",
    "        if int(p.x/10)<IB_MP_P_MYmap.shape[0]:\n",
    "            region = IB_MP_P_MYmap[int(p.x/10),int(p.y/10),int(p.z/10)]\n",
    "        else:\n",
    "            region=0\n",
    "        length = np.linalg.norm(np.array(edge.data[j].xyz)-np.array(edge.data[j+1].xyz))\n",
    "        lenMap[region]+=length\n",
    "    for k,v in lenMap.items():\n",
    "        clusterLength[edge.maxDepth][k]+=v\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.1267890713494533, 0.0, 0.0002458581360398717, 0.6297087866602614, 1.0]\n",
      "0.35134874322915094 0.39884307395603585 [6.78574433e-001 2.66540469e-001 1.94902976e-004 2.65315994e-010\n",
      " 6.72354074e-019 3.17192053e-030 2.78571192e-044 4.55449002e-061\n",
      " 1.38622254e-080 7.85445045e-103 8.28491966e-128 1.62686099e-155\n",
      " 5.94705888e-186 4.04709761e-219 5.12713518e-255]\n",
      "[0.2360189345362709, 1.0, 0.009141924587550377, 0.0012882537314322137, 0.0]\n",
      "0.2492898225710507 0.3860215262533808 [8.38953168e-001 1.55972348e-001 3.53082554e-005 9.73248336e-012\n",
      " 3.26655398e-021 1.33498003e-033 6.64321726e-049 4.02532456e-067\n",
      " 2.96990194e-088 2.66809984e-112 2.91864183e-139 3.88757058e-169\n",
      " 6.30513706e-202 1.24517254e-237 2.99421719e-276]\n",
      "[0.01678437860977126, 0.00014036518817409335, 0.03530458172980044, 1.0001403651881742, 0.0002448980186736376]\n",
      "0.2105229177469187 0.39502138274508963 [8.76221845e-001 1.37070469e-001 3.53249978e-005 1.49978433e-011\n",
      " 1.04902032e-020 1.20878041e-032 2.29466528e-047 7.17628537e-065\n",
      " 3.69733224e-085 3.13823656e-108 4.38824938e-134 1.01089336e-162\n",
      " 3.83643471e-194 2.39860481e-228 2.47057482e-265]\n",
      "[0.29563740979744296, 0.3376406074280959, 1.0, 0.0203253589659441, 0.0]\n",
      "0.33072067523829657 0.3619212491476219 [7.26063618e-001 1.99397579e-001 2.64795603e-005 1.70038530e-012\n",
      " 5.27994745e-023 7.92788866e-037 5.75613763e-054 2.02092632e-074\n",
      " 3.43095859e-098 2.81660637e-125 1.11810483e-155 2.14627107e-189\n",
      " 1.99219711e-226 8.94181591e-267 1.94072870e-310]\n",
      "[1.0, 0.45712954905055725, 0.0006289746788104211, 0.0, 0.0]\n",
      "0.2915517047458735 0.3959686445978235 [7.68290692e-001 2.03304253e-001 9.13852519e-005 6.97773575e-011\n",
      " 9.05027303e-020 1.99396532e-031 7.46246446e-046 4.74412046e-063\n",
      " 5.12315973e-083 9.39785722e-106 2.92838750e-131 1.55001894e-159\n",
      " 1.39365161e-190 2.12853204e-224 5.52223546e-261]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2oAAANOCAYAAABgFv8rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABFGUlEQVR4nO39f6xkh10e/j9P7bj8ComEtxXyj6wB02ICTcjKBCVVrS8g2UGyqZoW+ytUgtK4FbgChSK5pUqDq1akaSmf9GtKXUhDaRvHRLTakkX+VG0oKG2C1yWE2MZ049J6XaQsSYhEkuKavr9/3DG9LLt7J/bszsmd10u68pxz3jPzrHx0dp85Z87tzAQAAIDl+CPbDgAAAMAfpKgBAAAsjKIGAACwMIoaAADAwihqAAAAC3P5tt74yiuvnKNHj27r7QEAALbq4Ycf/q2ZOXKubVsrakePHs3Jkye39fYAAABb1fa/n2+bSx8BAAAWRlEDAABYmAOLWtu3t/1o2w+fZ3vbvq3tqbYfavt1m48JAACwO9Y5o/aOJDdfYPstSa5f/dyZ5B8//1gAAAC768CiNjO/kOTjFxi5Lck/nz3vT/Litl+6qYAAAAC7ZhN3fbwqyZP7lk+v1v3m2YNt78zeWbdce+21G3hrAABYz71/5T9sOwJb8t0/9v/ZdoTP2iW9mcjM3Dczx2bm2JEj5/x1AQAAADtvE0XtqSTX7Fu+erUOAACA52ATRe14kr+4uvvjK5N8cmb+0GWPAAAArOfA76i1fWeSm5Jc2fZ0kr+V5AVJMjM/luREktckOZXk00m+82KFBQAA2AUHFrWZueOA7ZPkuzeWCAAAYMdd0puJAAAAcDBFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBh1ipqbW9u+3jbU23vPsf2a9u+t+0vt/1Q29dsPioAAMBuOLCotb0syb1JbklyQ5I72t5w1tjfTPLAzLw8ye1JfnTTQQEAAHbFOmfUbkxyamaemJmnk9yf5LazZibJF68evyjJ/9xcRAAAgN2yTlG7KsmT+5ZPr9bt9+Yk3972dJITSf7quV6o7Z1tT7Y9eebMmecQFwAA4PDb1M1E7kjyjpm5OslrkvxU2z/02jNz38wcm5ljR44c2dBbAwAAHC7rFLWnklyzb/nq1br9Xp/kgSSZmf+c5POSXLmJgAAAALtmnaL2UJLr217X9ors3Szk+Fkz/yPJNyZJ26/KXlFzbSMAAMBzcGBRm5lnktyV5MEkj2Xv7o6PtL2n7a2rse9L8oa2v5LknUleNzNzsUIDAAAcZpevMzQzJ7J3k5D969607/GjSV612WgAAAC7aVM3EwEAAGBDFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYS5fZ6jtzUn+nySXJfnxmfmhc8z8hSRvTjJJfmVm/r8bzAkAHCKP/cmv2nYEtuSrfu2xbUeAzwkHFrW2lyW5N8k3Jzmd5KG2x2fm0X0z1yf560leNTOfaPvHLlZgAACAw26dSx9vTHJqZp6YmaeT3J/ktrNm3pDk3pn5RJLMzEc3GxMAAGB3rFPUrkry5L7l06t1+31lkq9s+762719dKgkAAMBzsNZ31NZ8neuT3JTk6iS/0PZrZua39w+1vTPJnUly7bXXbuitAQAADpd1zqg9leSafctXr9btdzrJ8Zn53zPz35L8evaK2x8wM/fNzLGZOXbkyJHnmhkAAOBQW6eoPZTk+rbXtb0iye1Jjp8182+ydzYtba/M3qWQT2wuJgAAwO44sKjNzDNJ7kryYJLHkjwwM4+0vaftrauxB5N8rO2jSd6b5Ptn5mMXKzQAAMBhttZ31GbmRJITZ617077Hk+SNqx8AAACeh3UufQQAAOASUtQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFmatotb25raPtz3V9u4LzP25ttP22OYiAgAA7JYDi1rby5Lcm+SWJDckuaPtDeeYe2GS70nygU2HBAAA2CXrnFG7McmpmXliZp5Ocn+S284x97eTvCXJ/9pgPgAAgJ2zTlG7KsmT+5ZPr9b9vrZfl+SamXnPhV6o7Z1tT7Y9eebMmc86LAAAwC543jcTaftHkvxwku87aHZm7puZYzNz7MiRI8/3rQEAAA6ldYraU0mu2bd89Wrds16Y5KVJfr7tbyR5ZZLjbigCAADw3KxT1B5Kcn3b69pekeT2JMef3Tgzn5yZK2fm6MwcTfL+JLfOzMmLkhgAAOCQO7CozcwzSe5K8mCSx5I8MDOPtL2n7a0XOyAAAMCuuXydoZk5keTEWevedJ7Zm55/LAAAgN31vG8mAgAAwGYpagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALs1ZRa3tz28fbnmp79zm2v7Hto20/1Pbft33J5qMCAADshgOLWtvLktyb5JYkNyS5o+0NZ439cpJjM/O1Sd6d5O9tOigAAMCuWOeM2o1JTs3MEzPzdJL7k9y2f2Bm3jszn14tvj/J1ZuNCQAAsDvWKWpXJXly3/Lp1brzeX2SnzvXhrZ3tj3Z9uSZM2fWTwkAALBDNnozkbbfnuRYkreea/vM3Dczx2bm2JEjRzb51gAAAIfG5WvMPJXkmn3LV6/W/QFtvynJDyT5MzPzu5uJBwAAsHvWOaP2UJLr217X9ooktyc5vn+g7cuT/JMkt87MRzcfEwAAYHccWNRm5pkkdyV5MMljSR6YmUfa3tP21tXYW5N8UZKfbvvBtsfP83IAAAAcYJ1LHzMzJ5KcOGvdm/Y9/qYN5wIAANhZG72ZCAAAAM+fogYAALAwihoAAMDCKGoAAAALs9bNRHbJ0bvfs+0IbMlv/NC3bDsCAAAkcUYNAABgcRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYdYqam1vbvt421Nt7z7H9j/a9l2r7R9oe3TjSQEAAHbEgUWt7WVJ7k1yS5IbktzR9oazxl6f5BMz8xVJ/mGSt2w6KAAAwK5Y54zajUlOzcwTM/N0kvuT3HbWzG1JfnL1+N1JvrFtNxcTAABgd1y+xsxVSZ7ct3w6ydefb2Zmnmn7ySRfkuS39g+1vTPJnavF32n7+HMJzUV1Zc76/7Yr6jzwNu3sfsfW2ffYht3e73yWv007u+/d9U+2neC8XnK+DesUtY2ZmfuS3Hcp35PPTtuTM3Ns2znYLfY7tsW+xzbY79gW+97nlnUufXwqyTX7lq9erTvnTNvLk7woycc2ERAAAGDXrFPUHkpyfdvr2l6R5PYkx8+aOZ7kO1aPX5vkP8zMbC4mAADA7jjw0sfVd87uSvJgksuSvH1mHml7T5KTM3M8yU8k+am2p5J8PHtljs9NLk1lG+x3bIt9j22w37Et9r3PIXXiCwAAYFnW+oXXAAAAXDqKGgAAwMIoajum7Yvbfte+5Zva/uw2M7E72k7bf7Fv+fK2Z9r+bNuvbvvrbT9/3/b3tL1jO2k5zNr+XtsPtv1w259u+wXbzsThc6Fj3mr5davlD7Z9tO0btpeWw6jt76z+e7TtZ1b72q+0/U9t/8S283FhitrueXGS7zpoaF2rX8cA6/pUkpfuK2PfnNWv+5iZR5L8TJIfSJK235rkBTPzzi3k5PD7zMy8bGZemuTpJH9l24E4lM57zNvnXTPzsiQ3Jfm7bf/4pYvHjvnI6rj3p5L8ZJK/se1AXJiidsi1fePqE+MPt/3eJD+U5MtXn6i8dTX2RW3f3fbX2v7Ltl099xVt/2Pbh9s+2PZLV+t/vu2PtD2Z5Hu28gfjc9mJJN+yenxHkv1F7J4kf77ty7K3r373pY3GjvrFJF+x7RAcWhc65v2+mfloko8kecklysVu++Ikn9h2CC5MUTvE2r4iyXcm+fokr0zyhiRvyf/9ROX7V6MvT/K9SW5I8mVJXtX2BUn+UZLXzswrkrw9yd/Z9/JXzMyxmfkHl+QPw2Fyf5Lb235ekq9N8oFnN8zMp5P8tSS/kOT+mfmv24nIrlhdFXBLkl/ddhYOrfMe8/Zr+2XZ+zv41CXMxm559oP6jyR5Y5If3nYgLsxla4fbq5P865n5VJK0/Zkkf/occ780M6dXMx9McjTJbyd5aZJ/tzrBdlmS39z3nHddrNAcbjPzobZHs/fJ8olzbP+3bX87yY9e4mjsls9fHe+SvTNqP7HFLBxiBx3zknxb21cn+d0kf3lmPn4p87FTPrK6zDZtvy17v1Pt5q0m4oIUNZK9vxye9XvZ2y+a5JGZ+YbzPOdTFz0Vh9nxJH8/e9/J+JJzbP8/qx+4WD7z7D9Y4BK40DHvXTNz1yVPxK47nuSfbTsEF+bSx8PtF5N8a9svaPuFSf5skvcleeEaz308yZG235AkbV/Q9qsvXlR2zNuT/ODMuNwM2AWOeSzNq7P3nUgWzBm1Q2xm/kvbdyT5pdWqH5+Zh9u+r+2Hk/xckvec57lPt31tkre1fVH29pUfSfLIxU/OYbe61PZt284BcCk45rEQX7665LvZu9vtX9puHA7Smdl2BgAAAPZx6SMAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACzM5dt64yuvvHKOHj26rbcHAADYqocffvi3ZubIubZtragdPXo0J0+e3NbbAwAAbFXb/36+bS59BAAAWBhFDQAAYGEOLGpt3972o20/fJ7tbfu2tqfafqjt120+JgAAwO5Y54zaO5LcfIHttyS5fvVzZ5J//PxjAQAA7K4Di9rM/EKSj19g5LYk/3z2vD/Ji9t+6aYCAgAA7JpN3PXxqiRP7ls+vVr3m2cPtr0ze2fdcu21127greEQefOLtp2AbXjzJ7edAABYoEt6M5GZuW9mjs3MsSNHzvnrAgAAAHbeJoraU0mu2bd89WodAAAAz8EmitrxJH9xdffHVyb55Mz8ocseAQAAWM+B31Fr+84kNyW5su3pJH8ryQuSZGZ+LMmJJK9JcirJp5N858UKCwAAsAsOLGozc8cB2yfJd28sEQAAwI67pDcTAQAA4GCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwaxW1tje3fbztqbZ3n2P7tW3f2/aX236o7Ws2HxUAAGA3HFjU2l6W5N4ktyS5IckdbW84a+xvJnlgZl6e5PYkP7rpoAAAALtinTNqNyY5NTNPzMzTSe5PcttZM5Pki1ePX5Tkf24uIgAAwG65fI2Zq5I8uW/5dJKvP2vmzUn+37Z/NckXJvmmjaQDAADYQZu6mcgdSd4xM1cneU2Sn2r7h1677Z1tT7Y9eebMmQ29NQAAwOGyTlF7Ksk1+5avXq3b7/VJHkiSmfnPST4vyZVnv9DM3Dczx2bm2JEjR55bYgAAgENunaL2UJLr217X9ors3Szk+Fkz/yPJNyZJ26/KXlFzygwAAOA5OLCozcwzSe5K8mCSx7J3d8dH2t7T9tbV2PcleUPbX0nyziSvm5m5WKEBAAAOs3VuJpKZOZHkxFnr3rTv8aNJXrXZaAAAALtpUzcTAQAAYEMUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFWauotb257eNtT7W9+zwzf6Hto20fafuvNhsTAABgd1x+0EDby5Lcm+Sbk5xO8lDb4zPz6L6Z65P89SSvmplPtP1jFyswAADAYbfOGbUbk5yamSdm5ukk9ye57ayZNyS5d2Y+kSQz89HNxgQAANgd6xS1q5I8uW/59Grdfl+Z5Cvbvq/t+9vefK4Xantn25NtT545c+a5JQYAADjkNnUzkcuTXJ/kpiR3JPmnbV989tDM3Dczx2bm2JEjRzb01gAAAIfLOkXtqSTX7Fu+erVuv9NJjs/M/56Z/5bk17NX3AAAAPgsrVPUHkpyfdvr2l6R5PYkx8+a+TfZO5uWtldm71LIJzYXEwAAYHccWNRm5pkkdyV5MMljSR6YmUfa3tP21tXYg0k+1vbRJO9N8v0z87GLFRoAAOAwO/D2/EkyMyeSnDhr3Zv2PZ4kb1z9AAAA8Dxs6mYiAAAAbIiiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwaxW1tje3fbztqbZ3X2Duz7Wdtsc2FxEAAGC3HFjU2l6W5N4ktyS5IckdbW84x9wLk3xPkg9sOiQAAMAuWeeM2o1JTs3MEzPzdJL7k9x2jrm/neQtSf7XBvMBAADsnHWK2lVJnty3fHq17ve1/bok18zMey70Qm3vbHuy7ckzZ8581mEBAAB2wfO+mUjbP5Lkh5N830GzM3PfzBybmWNHjhx5vm8NAABwKK1T1J5Kcs2+5atX6571wiQvTfLzbX8jySuTHHdDEQAAgOdmnaL2UJLr217X9ooktyc5/uzGmfnkzFw5M0dn5miS9ye5dWZOXpTEAAAAh9yBRW1mnklyV5IHkzyW5IGZeaTtPW1vvdgBAQAAds3l6wzNzIkkJ85a96bzzN70/GMBAADsrud9MxEAAAA2S1EDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhbl82wGW5ujd79l2BLbkN37oW7YdAQAAkjijBgAAsDhrFbW2N7d9vO2ptnefY/sb2z7a9kNt/33bl2w+KgAAwG44sKi1vSzJvUluSXJDkjva3nDW2C8nOTYzX5vk3Un+3qaDAgAA7Ip1zqjdmOTUzDwxM08nuT/JbfsHZua9M/Pp1eL7k1y92ZgAAAC7Y52idlWSJ/ctn16tO5/XJ/m5c21oe2fbk21PnjlzZv2UAAAAO2SjNxNp++1JjiV567m2z8x9M3NsZo4dOXJkk28NAABwaKxze/6nklyzb/nq1bo/oO03JfmBJH9mZn53M/EAAAB2zzpn1B5Kcn3b69pekeT2JMf3D7R9eZJ/kuTWmfno5mMCAADsjgOL2sw8k+SuJA8meSzJAzPzSNt72t66Gntrki9K8tNtP9j2+HleDgAAgAOsc+ljZuZEkhNnrXvTvsfftOFcAAAAO2ujNxMBAADg+VPUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVZq6i1vbnt421Ptb37HNv/aNt3rbZ/oO3RjScFAADYEQcWtbaXJbk3yS1JbkhyR9sbzhp7fZJPzMxXJPmHSd6y6aAAAAC74vI1Zm5McmpmnkiStvcnuS3Jo/tmbkvy5tXjdyf5/7XtzMwGswKwYV/zk1+z7Qhsya9+x69uOwIAF7DOpY9XJXly3/Lp1bpzzszMM0k+meRLNhEQAABg16xzRm1j2t6Z5M7V4u+0ffxSvj9ruTLJb207xDbUBbvbtLP7XX6w206w63Z23+vr7HtbtLP7HVtn31uel5xvwzpF7akk1+xbvnq17lwzp9tenuRFST529gvNzH1J7lvjPdmStidn5ti2c7Bb7Hdsi32PbbDfsS32vc8t61z6+FCS69te1/aKJLcnOX7WzPEk37F6/Nok/8H30wAAAJ6bA8+ozcwzbe9K8mCSy5K8fWYeaXtPkpMzczzJTyT5qbanknw8e2UOAACA52Ct76jNzIkkJ85a96Z9j/9Xkj+/2WhsiUtT2Qb7Hdti32Mb7Hdsi33vc0hdoQgAALAs63xHDQAAgEtIUdsxbV/c9rv2Ld/U9me3mYnd0Xba/ot9y5e3PdP2Z9t+ddtfb/v5+7a/p+0d20nLYdb299p+sO2H2/502y/YdiYOnwsd81bLr1stf7Dto23fsL20HEZtf2f136NtP7Pa136l7X9q+ye2nY8LU9R2z4uTfNdBQ+ta/ToGWNenkrx0Xxn75qx+3cfMPJLkZ5L8QJK0/dYkL5iZd24hJ4ffZ2bmZTPz0iRPJ/kr2w7EoXTeY94+75qZlyW5KcnfbfvHL108dsxHVse9P5XkJ5P8jW0H4sIUtUOu7RtXnxh/uO33JvmhJF+++kTlrauxL2r77ra/1vZftu3qua9o+x/bPtz2wbZfulr/821/pO3JJN+zlT8Yn8tOJPmW1eM7kuwvYvck+fNtX5a9ffW7L200dtQvJvmKbYfg0LrQMe/3zcxHk3wkF/jlt7BBX5zkE9sOwYUpaodY21ck+c4kX5/klUnekOQt+b+fqHz/avTlSb43yQ1JvizJq9q+IMk/SvLamXlFkrcn+Tv7Xv6KmTk2M//gkvxhOEzuT3J7289L8rVJPvDshpn5dJK/luQXktw/M/91OxHZFaurAm5J8qvbzsKhdd5j3n5tvyx7fwefuoTZ2C3PflD/kSRvTPLD2w7Ehbls7XB7dZJ/PTOfSpK2P5PkT59j7pdm5vRq5oNJjib57SQvTfLvVifYLkvym/ue866LFZrDbWY+1PZo9j5ZPnGO7f+27W8n+dFLHI3d8vmr412yd0btJ7aYhUPsoGNekm9r++okv5vkL8/Mxy9lPnbKR1aX2abtt2XvVv03bzURF6Sokez95fCs38veftEkj8zMN5znOZ+66Kk4zI4n+fvZ+07Gl5xj+/9Z/cDF8pln/8ECl8CFjnnvmpm7Lnkidt3xJP9s2yG4MJc+Hm6/mORb235B2y9M8meTvC/JC9d47uNJjrT9hiRp+4K2X33xorJj3p7kB2fG5WbALnDMY2lenb3vRLJgzqgdYjPzX9q+I8kvrVb9+Mw83PZ9bT+c5OeSvOc8z3267WuTvK3ti7K3r/xIkkcufnIOu9Wltm/bdg6AS8Exj4X48tUl383e3W7/0nbjcJDOzLYzAAAAsI9LHwEAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIW5fFtvfOWVV87Ro0e39fYAAABb9fDDD//WzBw517atFbWjR4/m5MmT23p7AACArWr738+3zaWPAAAAC3NgUWv79rYfbfvh82xv27e1PdX2Q22/bvMxAQAAdsc6Z9TekeTmC2y/Jcn1q587k/zj5x8LAABgdx1Y1GbmF5J8/AIjtyX557Pn/Ule3PZLNxUQAABg12ziO2pXJXly3/Lp1ToAAACeg0t618e2d2bv8shce+21l/KtAYAFeexPftW2I7AlX/Vrj207AnxO2MQZtaeSXLNv+erVuj9kZu6bmWMzc+zIkXP+ugAAAICdt4midjzJX1zd/fGVST45M7+5gdcFAADYSQde+tj2nUluSnJl29NJ/laSFyTJzPxYkhNJXpPkVJJPJ/nOixUWAABgFxxY1GbmjgO2T5Lv3lgiAACAHbeJSx8BAADYIEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGHWKmptb277eNtTbe8+x/Zr27637S+3/VDb12w+KgAAwG44sKi1vSzJvUluSXJDkjva3nDW2N9M8sDMvDzJ7Ul+dNNBAQAAdsU6Z9RuTHJqZp6YmaeT3J/ktrNmJskXrx6/KMn/3FxEAACA3bJOUbsqyZP7lk+v1u335iTf3vZ0khNJ/uq5XqjtnW1Ptj155syZ5xAXAADg8NvUzUTuSPKOmbk6yWuS/FTbP/TaM3PfzBybmWNHjhzZ0FsDAAAcLusUtaeSXLNv+erVuv1en+SBJJmZ/5zk85JcuYmAAAAAu2adovZQkuvbXtf2iuzdLOT4WTP/I8k3Jknbr8peUXNtIwAAwHNwYFGbmWeS3JXkwSSPZe/ujo+0vaftraux70vyhra/kuSdSV43M3OxQgMAABxml68zNDMnsneTkP3r3rTv8aNJXrXZaAAAALtpUzcTAQAAYEMUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFWauotb257eNtT7W9+zwzf6Hto20fafuvNhsTAABgd1x+0EDby5Lcm+Sbk5xO8lDb4zPz6L6Z65P89SSvmplPtP1jFyswAADAYbfOGbUbk5yamSdm5ukk9ye57ayZNyS5d2Y+kSQz89HNxgQAANgd6xS1q5I8uW/59Grdfl+Z5Cvbvq/t+9vefK4Xantn25NtT545c+a5JQYAADjkNnUzkcuTXJ/kpiR3JPmnbV989tDM3Dczx2bm2JEjRzb01gAAAIfLOkXtqSTX7Fu+erVuv9NJjs/M/56Z/5bk17NX3AAAAPgsrVPUHkpyfdvr2l6R5PYkx8+a+TfZO5uWtldm71LIJzYXEwAAYHccWNRm5pkkdyV5MMljSR6YmUfa3tP21tXYg0k+1vbRJO9N8v0z87GLFRoAAOAwO/D2/EkyMyeSnDhr3Zv2PZ4kb1z9AAAA8Dxs6mYiAAAAbIiiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwaxW1tje3fbztqbZ3X2Duz7Wdtsc2FxEAAGC3HFjU2l6W5N4ktyS5IckdbW84x9wLk3xPkg9sOiQAAMAuWeeM2o1JTs3MEzPzdJL7k9x2jrm/neQtSf7XBvMBAADsnHWK2lVJnty3fHq17ve1/bok18zMey70Qm3vbHuy7ckzZ8581mEBAAB2wfO+mUjbP5Lkh5N830GzM3PfzBybmWNHjhx5vm8NAABwKK1T1J5Kcs2+5atX6571wiQvTfLzbX8jySuTHHdDEQAAgOdmnaL2UJLr217X9ooktyc5/uzGmfnkzFw5M0dn5miS9ye5dWZOXpTEAAAAh9yBRW1mnklyV5IHkzyW5IGZeaTtPW1vvdgBAQAAds3l6wzNzIkkJ85a96bzzN70/GMBAADsrud9MxEAAAA2S1EDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFmatotb25raPtz3V9u5zbH9j20fbfqjtv2/7ks1HBQAA2A0HFrW2lyW5N8ktSW5IckfbG84a++Ukx2bma5O8O8nf23RQAACAXbHOGbUbk5yamSdm5ukk9ye5bf/AzLx3Zj69Wnx/kqs3GxMAAGB3rFPUrkry5L7l06t15/P6JD/3fEIBAADssss3+WJtvz3JsSR/5jzb70xyZ5Jce+21m3xrAACAQ2OdM2pPJblm3/LVq3V/QNtvSvIDSW6dmd891wvNzH0zc2xmjh05cuS55AUAADj01ilqDyW5vu11ba9IcnuS4/sH2r48yT/JXkn76OZjAgAA7I4Di9rMPJPkriQPJnksyQMz80jbe9reuhp7a5IvSvLTbT/Y9vh5Xg4AAIADrPUdtZk5keTEWevetO/xN204FwAAwM5a6xdeAwAAcOkoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCXL7tAABsz9f85NdsOwJb8qvf8avbjgDABax1Rq3tzW0fb3uq7d3n2P5H275rtf0DbY9uPCkAAMCOOPCMWtvLktyb5JuTnE7yUNvjM/PovrHXJ/nEzHxF29uTvCXJt12MwBfb0bvfs+0IbMlv/NC3bDsCAAAkWe+M2o1JTs3MEzPzdJL7k9x21sxtSX5y9fjdSb6xbTcXEwAAYHes8x21q5I8uW/5dJKvP9/MzDzT9pNJviTJb+0fantnkjtXi7/T9vHnEpqL6sqc9f9tV/Qt206w03Z2v2Prdnbf6+t8nrpFO7vfJUl8lr9Nu73vLdNLzrfhkt5MZGbuS3LfpXxPPjttT87MsW3nYLfY79gW+x7bYL9jW+x7n1vWufTxqSTX7Fu+erXunDNtL0/yoiQf20RAAACAXbNOUXsoyfVtr2t7RZLbkxw/a+Z4ku9YPX5tkv8wM7O5mAAAALvjwEsfV985uyvJg0kuS/L2mXmk7T1JTs7M8SQ/keSn2p5K8vHslTk+N7k0lW2w37Et9j22wX7Httj3PofUiS8AAIBlWesXXgMAAHDpKGoAAAALo6jtmLYvbvtd+5Zvavuz28zE7mg7bf/FvuXL255p+7Ntv7rtr7f9/H3b39P2ju2k5TBr+3ttP9j2w21/uu0XbDsTh8+Fjnmr5detlj/Y9tG2b9heWg6jtr+z+u/Rtp9Z7Wu/0vY/tf0T287HhSlqu+fFSb7roKF1rX4dA6zrU0leuq+MfXNWv+5jZh5J8jNJfiBJ2n5rkhfMzDu3kJPD7zMz87KZeWmSp5P8lW0H4lA67zFvn3fNzMuS3JTk77b945cuHjvmI6vj3p9K8pNJ/sa2A3Fhitoh1/aNq0+MP9z2e5P8UJIvX32i8tbV2Be1fXfbX2v7L9t29dxXtP2PbR9u+2DbL12t//m2P9L2ZJLv2cofjM9lJ5J8y+rxHUn2F7F7kvz5ti/L3r763Zc2GjvqF5N8xbZDcGhd6Jj3+2bmo0k+kuQllygXu+2Lk3xi2yG4MEXtEGv7iiTfmeTrk7wyyRuSvCX/9xOV71+NvjzJ9ya5IcmXJXlV2xck+UdJXjszr0jy9iR/Z9/LXzEzx2bmH1ySPwyHyf1Jbm/7eUm+NskHnt0wM59O8teS/EKS+2fmv24nIrtidVXALUl+ddtZOLTOe8zbr+2XZe/v4FOXMBu75dkP6j+S5I1Jfnjbgbgwl60dbq9O8q9n5lNJ0vZnkvzpc8z90sycXs18MMnRJL+d5KVJ/t3qBNtlSX5z33PedbFCc7jNzIfaHs3eJ8snzrH937b97SQ/eomjsVs+f3W8S/bOqP3EFrNwiB10zEvybW1fneR3k/zlmfn4pczHTvnI6jLbtP227P1OtZu3mogLUtRI9v5yeNbvZW+/aJJHZuYbzvOcT130VBxmx5P8/ex9J+NLzrH9/6x+4GL5zLP/YIFL4ELHvHfNzF2XPBG77niSf7btEFyYSx8Pt19M8q1tv6DtFyb5s0nel+SFazz38SRH2n5DkrR9QduvvnhR2TFvT/KDM+NyM2AXOOaxNK/O3nciWTBn1A6xmfkvbd+R5JdWq358Zh5u+762H07yc0nec57nPt32tUne1vZF2dtXfiTJIxc/OYfd6lLbt207B8Cl4JjHQnz56pLvZu9ut39pu3E4SGdm2xkAAADYx6WPAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwly+rTe+8sor5+jRo9t6ewAAgK16+OGHf2tmjpxr29aK2tGjR3Py5MltvT0AAMBWtf3v59vm0kcAAICFObCotX1724+2/fB5trft29qeavuhtl+3+ZgAAAC7Y50zau9IcvMFtt+S5PrVz51J/vHzjwUAALC7DixqM/MLST5+gZHbkvzz2fP+JC9u+6WbCggAALBrNnEzkauSPLlv+fRq3W+ePdj2zuyddcu11167gbcG4Pn4mp/8mm1HYEt+9Tt+ddsRALiAS3ozkZm5b2aOzcyxI0fOeRdKAACAnbeJovZUkmv2LV+9WgcAAMBzsImidjzJX1zd/fGVST45M3/oskcAAADWc+B31Nq+M8lNSa5sezrJ30rygiSZmR9LciLJa5KcSvLpJN95scICAADsggOL2szcccD2SfLdG0sEAACw4y7pzUQAAAA4mKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDBrFbW2N7d9vO2ptnefY/u1bd/b9pfbfqjtazYfFQAAYDccWNTaXpbk3iS3JLkhyR1tbzhr7G8meWBmXp7k9iQ/uumgAAAAu2KdM2o3Jjk1M0/MzNNJ7k9y21kzk+SLV49flOR/bi4iAADAblmnqF2V5Ml9y6dX6/Z7c5Jvb3s6yYkkf/VcL9T2zrYn2548c+bMc4gLAABw+G3qZiJ3JHnHzFyd5DVJfqrtH3rtmblvZo7NzLEjR45s6K0BAAAOl3WK2lNJrtm3fPVq3X6vT/JAkszMf07yeUmu3ERAAACAXbNOUXsoyfVtr2t7RfZuFnL8rJn/keQbk6TtV2WvqLm2EQAA4Dk4sKjNzDNJ7kryYJLHsnd3x0fa3tP21tXY9yV5Q9tfSfLOJK+bmblYoQEAAA6zy9cZmpkT2btJyP51b9r3+NEkr9psNAAAgN20qZuJAAAAsCGKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCrFXU2t7c9vG2p9refZ6Zv9D20baPtP1Xm40JAACwOy4/aKDtZUnuTfLNSU4neajt8Zl5dN/M9Un+epJXzcwn2v6xixUYAADgsFvnjNqNSU7NzBMz83SS+5PcdtbMG5LcOzOfSJKZ+ehmYwIAAOyOdYraVUme3Ld8erVuv69M8pVt39f2/W1vPtcLtb2z7cm2J8+cOfPcEgMAABxym7qZyOVJrk9yU5I7kvzTti8+e2hm7puZYzNz7MiRIxt6awAAgMNlnaL2VJJr9i1fvVq33+kkx2fmf8/Mf0vy69krbgAAAHyW1ilqDyW5vu11ba9IcnuS42fN/JvsnU1L2yuzdynkE5uLCQAAsDsOLGoz80ySu5I8mOSxJA/MzCNt72l762rswSQfa/tokvcm+f6Z+djFCg0AAHCYHXh7/iSZmRNJTpy17k37Hk+SN65+AAAAeB42dTMRAAAANkRRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZmraLW9ua2j7c91fbuC8z9ubbT9tjmIgIAAOyWA4ta28uS3JvkliQ3JLmj7Q3nmHthku9J8oFNhwQAANgl65xRuzHJqZl5YmaeTnJ/ktvOMfe3k7wlyf/aYD4AAICds05RuyrJk/uWT6/W/b62X5fkmpl5zwazAQAA7KTnfTORtn8kyQ8n+b41Zu9se7LtyTNnzjzftwYAADiU1ilqTyW5Zt/y1at1z3phkpcm+fm2v5HklUmOn+uGIjNz38wcm5ljR44cee6pAQAADrF1itpDSa5ve13bK5LcnuT4sxtn5pMzc+XMHJ2Zo0nen+TWmTl5URIDAAAccgcWtZl5JsldSR5M8liSB2bmkbb3tL31YgcEAADYNZevMzQzJ5KcOGvdm84ze9PzjwUAALC7nvfNRAAAANgsRQ0AAGBhFDUAAICFWes7asAl8OYXbTsB2/DmT247AQCwQM6oAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMuz6e5ejd79l2BLbkN37oW7YdAQAAkjijBgAAsDiKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAszFpFre3NbR9ve6rt3efY/sa2j7b9UNt/3/Ylm48KAACwGw4sam0vS3JvkluS3JDkjrY3nDX2y0mOzczXJnl3kr+36aAAAAC7Yp0zajcmOTUzT8zM00nuT3Lb/oGZee/MfHq1+P4kV282JgAAwO5Yp6hdleTJfcunV+vO5/VJfu5cG9re2fZk25NnzpxZPyUAAMAO2ejNRNp+e5JjSd56ru0zc9/MHJuZY0eOHNnkWwMAABwal68x81SSa/YtX71a9we0/aYkP5Dkz8zM724mHgAAwO5Z54zaQ0mub3td2yuS3J7k+P6Bti9P8k+S3DozH918TAAAgN1xYFGbmWeS3JXkwSSPJXlgZh5pe0/bW1djb03yRUl+uu0H2x4/z8sBAABwgHUufczMnEhy4qx1b9r3+Js2nAsAAGBnbfRmIgAAADx/ihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwly+zlDbm5P8P0kuS/LjM/NDZ23/o0n+eZJXJPlYkm+bmd/YbFQA4LB47E9+1bYjsCVf9WuPbTsCfE448Ixa28uS3JvkliQ3JLmj7Q1njb0+ySdm5iuS/MMkb9l0UAAAgF2xzqWPNyY5NTNPzMzTSe5PcttZM7cl+cnV43cn+ca23VxMAACA3bHOpY9XJXly3/LpJF9/vpmZeabtJ5N8SZLf2j/U9s4kd64Wf6ft488lNBfVlTnr/9uuqPPA27Sz+11+0GdaW7az+15fZ9/bop3d75IkPsvfpt3e95bpJefbsNZ31DZlZu5Lct+lfE8+O21Pzsyxbedgt9jv2Bb7Httgv2Nb7HufW9a59PGpJNfsW756te6cM20vT/Ki7N1UBAAAgM/SOkXtoSTXt72u7RVJbk9y/KyZ40m+Y/X4tUn+w8zM5mICAADsjgMvfVx95+yuJA9m7/b8b5+ZR9rek+TkzBxP8hNJfqrtqSQfz16Z43OTS1PZBvsd22LfYxvsd2yLfe9zSJ34AgAAWJZ1Ln0EAADgElLUAAAAFkZR2zFtX9z2u/Yt39T2Z7eZid3Rdtr+i33Ll7c90/Zn2351219v+/n7tr+n7R3bScth1vb32n6w7Yfb/nTbL9h2Jg6fCx3zVsuvWy1/sO2jbd+wvbQcRm1/Z/Xfo20/s9rXfqXtf2r7J7adjwtT1HbPi5N810FD61r9OgZY16eSvHRfGfvmrH7dx8w8kuRnkvxAkrT91iQvmJl3biEnh99nZuZlM/PSJE8n+SvbDsShdN5j3j7vmpmXJbkpyd9t+8cvXTx2zEdWx70/leQnk/yNbQfiwhS1Q67tG1efGH+47fcm+aEkX776ROWtq7Evavvutr/W9l+27eq5r2j7H9s+3PbBtl+6Wv/zbX+k7ckk37OVPxify04k+ZbV4zuS7C9i9yT5821flr199bsvbTR21C8m+Ypth+DQutAx7/fNzEeTfCTJSy5RLnbbFyf5xLZDcGGK2iHW9hVJvjPJ1yd5ZZI3JHlL/u8nKt+/Gn15ku9NckOSL0vyqrYvSPKPkrx2Zl6R5O1J/s6+l79iZo7NzD+4JH8YDpP7k9ze9vOSfG2SDzy7YWY+neSvJfmFJPfPzH/dTkR2xeqqgFuS/Oq2s3BonfeYt1/bL8ve38GnLmE2dsuzH9R/JMkbk/zwtgNxYS5bO9xeneRfz8ynkqTtzyT50+eY+6WZOb2a+WCSo0l+O8lLk/y71Qm2y5L85r7nvOtiheZwm5kPtT2avU+WT5xj+79t+9tJfvQSR2O3fP7qeJfsnVH7iS1m4RA76JiX5NvavjrJ7yb5yzPz8UuZj53ykdVltmn7bdn7nWo3bzURF6Sokez95fCs38veftEkj8zMN5znOZ+66Kk4zI4n+fvZ+07Gl5xj+/9Z/cDF8pln/8ECl8CFjnnvmpm7Lnkidt3xJP9s2yG4MJc+Hm6/mORb235B2y9M8meTvC/JC9d47uNJjrT9hiRp+4K2X33xorJj3p7kB2fG5WbALnDMY2lenb3vRLJgzqgdYjPzX9q+I8kvrVb9+Mw83PZ9bT+c5OeSvOc8z3267WuTvK3ti7K3r/xIkkcufnIOu9Wltm/bdg6AS8Exj4X48tUl383e3W7/0nbjcJDOzLYzAAAAsI9LHwEAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEu39YbX3nllXP06NFtvT0AAMBWPfzww781M0fOtW1rRe3o0aM5efLktt4eAABgq9r+9/Ntc+kjAADAwihqAAAAC3NgUWv79rYfbfvh82xv27e1PdX2Q22/bvMxAQAAdsc6Z9TekeTmC2y/Jcn1q587k/zj5x8LAABgdx1Y1GbmF5J8/AIjtyX557Pn/Ule3PZLNxUQAABg12ziro9XJXly3/Lp1brfPHuw7Z3ZO+uWa6+9dgNvvXlH737PtiOwJb/xQ9+y7QgAAJDkEt9MZGbum5ljM3PsyJFz/roAAACAnbeJovZUkmv2LV+9WgcAAMBzsImidjzJX1zd/fGVST45M3/oskcAAADWc+B31Nq+M8lNSa5sezrJ30rygiSZmR9LciLJa5KcSvLpJN95scICAADsggOL2szcccD2SfLdG0sEAACw4y7pzUQAAAA4mKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALMxaRa3tzW0fb3uq7d3n2H5t2/e2/eW2H2r7ms1HBQAA2A0HFrW2lyW5N8ktSW5IckfbG84a+5tJHpiZlye5PcmPbjooAADArljnjNqNSU7NzBMz83SS+5PcdtbMJPni1eMXJfmfm4sIAACwWy5fY+aqJE/uWz6d5OvPmnlzkv+37V9N8oVJvmkj6QAAAHbQpm4mckeSd8zM1Ulek+Sn2v6h1257Z9uTbU+eOXNmQ28NAABwuKxT1J5Kcs2+5atX6/Z7fZIHkmRm/nOSz0ty5dkvNDP3zcyxmTl25MiR55YYAADgkFunqD2U5Pq217W9Ins3Czl+1sz/SPKNSdL2q7JX1JwyAwAAeA4OLGoz80ySu5I8mOSx7N3d8ZG297S9dTX2fUne0PZXkrwzyetmZi5WaAAAgMNsnZuJZGZOJDlx1ro37Xv8aJJXbTYaAADAbtrUzUQAAADYEEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGHWKmptb277eNtTbe8+z8xfaPto20fa/qvNxgQAANgdlx800PayJPcm+eYkp5M81Pb4zDy6b+b6JH89yatm5hNt/9jFCgwAAHDYrXNG7cYkp2bmiZl5Osn9SW47a+YNSe6dmU8kycx8dLMxAQAAdsc6Re2qJE/uWz69WrffVyb5yrbva/v+tjef64Xa3tn2ZNuTZ86ceW6JAQAADrlN3Uzk8iTXJ7kpyR1J/mnbF589NDP3zcyxmTl25MiRDb01AADA4bJOUXsqyTX7lq9erdvvdJLjM/O/Z+a/Jfn17BU3AAAAPkvrFLWHklzf9rq2VyS5Pcnxs2b+TfbOpqXtldm7FPKJzcUEAADYHQcWtZl5JsldSR5M8liSB2bmkbb3tL11NfZgko+1fTTJe5N8/8x87GKFBgAAOMwOvD1/kszMiSQnzlr3pn2PJ8kbVz8AAAA8D5u6mQgAAAAboqgBAAAsjKIGAACwMIoaAADAwihqAAAAC6OoAQAALIyiBgAAsDCKGgAAwMIoagAAAAujqAEAACyMogYAALAwihoAAMDCKGoAAAALo6gBAAAszOXbDgCsvPlF207ANrz5k9tOAAAskDNqAAAAC6OoAQAALIyiBgAAsDBrFbW2N7d9vO2ptndfYO7PtZ22xzYXEQAAYLccWNTaXpbk3iS3JLkhyR1tbzjH3AuTfE+SD2w6JAAAwC5Z54zajUlOzcwTM/N0kvuT3HaOub+d5C1J/tcG8wEAAOycdYraVUme3Ld8erXu97X9uiTXzMx7LvRCbe9se7LtyTNnznzWYQEAAHbB876ZSNs/kuSHk3zfQbMzc9/MHJuZY0eOHHm+bw0AAHAorVPUnkpyzb7lq1frnvXCJC9N8vNtfyPJK5Mcd0MRAACA52adovZQkuvbXtf2iiS3Jzn+7MaZ+eTMXDkzR2fmaJL3J7l1Zk5elMQAAACH3IFFbWaeSXJXkgeTPJbkgZl5pO09bW+92AEBAAB2zeXrDM3MiSQnzlr3pvPM3vT8YwEAAOyu530zEQAAADZLUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYmLWKWtub2z7e9lTbu8+x/Y1tH237obb/vu1LNh8VAABgNxxY1NpeluTeJLckuSHJHW1vOGvsl5Mcm5mvTfLuJH9v00EBAAB2xTpn1G5McmpmnpiZp5Pcn+S2/QMz896Z+fRq8f1Jrt5sTAAAgN2xTlG7KsmT+5ZPr9adz+uT/Ny5NrS9s+3JtifPnDmzfkoAAIAdstGbibT99iTHkrz1XNtn5r6ZOTYzx44cObLJtwYAADg0Ll9j5qkk1+xbvnq17g9o+01JfiDJn5mZ391MPAAAgN2zzhm1h5Jc3/a6tlckuT3J8f0DbV+e5J8kuXVmPrr5mAAAALvjwKI2M88kuSvJg0keS/LAzDzS9p62t67G3prki5L8dNsPtj1+npcDAADgAOtc+piZOZHkxFnr3rTv8TdtOBcAAMDO2ujNRAAAAHj+FDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVHUAAAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBhFDUAAICFUdQAAAAWRlEDAABYGEUNAABgYRQ1AACAhVmrqLW9ue3jbU+1vfsc2/9o23ettn+g7dGNJwUAANgRBxa1tpcluTfJLUluSHJH2xvOGnt9kk/MzFck+YdJ3rLpoAAAALtinTNqNyY5NTNPzMzTSe5PcttZM7cl+cnV43cn+ca23VxMAACA3XH5GjNXJXly3/LpJF9/vpmZeabtJ5N8SZLf2j/U9s4kd64Wf6ft488lNBfVlTnr/9uuqPPA27Sz+11+0GdaW7a7+x7bZL9jW+x7y/OS821Yp6htzMzcl+S+S/mefHbanpyZY9vOwW6x37Et9j22wX7Httj3Presc+njU0mu2bd89WrdOWfaXp7kRUk+tomAAAAAu2adovZQkuvbXtf2iiS3Jzl+1szxJN+xevzaJP9hZmZzMQEAAHbHgZc+rr5zdleSB5NcluTtM/NI23uSnJyZ40l+IslPtT2V5OPZK3N8bnJpKttgv2Nb7Htsg/2ObbHvfQ6pE18AAADLstYvvAYAAODSUdQAAAAWRlHbMW1f3Pa79i3f1PZnt5mJ3dF22v6LfcuXtz3T9mfbfnXbX2/7+fu2v6ftHdtJy2HW9vfafrDth9v+dNsv2HYmDp8LHfNWy69bLX+w7aNt37C9tBxGbX9n9d+jbT+z2td+pe1/avsntp2PC1PUds+Lk3zXQUPrWv06BljXp5K8dF8Z++asft3HzDyS5GeS/ECStP3WJC+YmXduISeH32dm5mUz89IkTyf5K9sOxKF03mPePu+amZcluSnJ3237xy9dPHbMR1bHvT+V5CeT/I1tB+LCFLVDru0bV58Yf7jt9yb5oSRfvvpE5a2rsS9q++62v9b2X7bt6rmvaPsf2z7c9sG2X7pa//Ntf6TtySTfs5U/GJ/LTiT5ltXjO5LsL2L3JPnzbV+WvX31uy9tNHbULyb5im2H4NC60DHv983MR5N8JMlLLlEudtsXJ/nEtkNwYYraIdb2FUm+M8nXJ3llkjckeUv+7ycq378afXmS701yQ5IvS/Kqti9I8o+SvHZmXpHk7Un+zr6Xv2Jmjs3MP7gkfxgOk/uT3N7285J8bZIPPLthZj6d5K8l+YUk98/Mf91ORHbF6qqAW5L86razcGid95i3X9svy97fwacuYTZ2y7Mf1H8kyRuT/PC2A3FhLls73F6d5F/PzKeSpO3PJPnT55j7pZk5vZr5YJKjSX47yUuT/LvVCbbLkvzmvue862KF5nCbmQ+1PZq9T5ZPnGP7v23720l+9BJHY7d8/up4l+ydUfuJLWbhEDvomJfk29q+OsnvJvnLM/PxS5mPnfKR1WW2aftt2fudajdvNREXpKiR7P3l8Kzfy95+0SSPzMw3nOc5n7roqTjMjif5+9n7TsaXnGP7/1n9wMXymWf/wQKXwIWOee+ambsueSJ23fEk/2zbIbgwlz4ebr+Y5FvbfkHbL0zyZ5O8L8kL13ju40mOtP2GJGn7grZfffGismPenuQHZ8blZsAucMxjaV6dve9EsmDOqB1iM/Nf2r4jyS+tVv34zDzc9n1tP5zk55K85zzPfbrta5O8re2Lsrev/EiSRy5+cg671aW2b9t2DoBLwTGPhfjy1SXfzd7dbv/SduNwkM7MtjMAAACwj0sfAQAAFkZRAwAAWBhFDQAAYGEUNQAAgIVR1AAAABZGUQMAAFgYRQ0AAGBh/v8eY53HJDkzJAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1080x1080 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "clusterLength\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.mlab as mlab\n",
    "matplotlib.use('module://matplotlib_inline.backend_inline')\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "from scipy.stats import norm\n",
    "\n",
    "def Normalize(data):\n",
    "    m = np.mean(data)\n",
    "    mx = max(data)\n",
    "    mn = min(data)\n",
    "    return [(float(i)) / (mx - mn) for i in data]\n",
    "\n",
    "plt.figure(figsize=(15,15))\n",
    "for j in range(5):\n",
    "    plt.subplot(5,1,j+1)\n",
    "    x= Normalize(np.array(list(clusterLength[j].values())))\n",
    "    print(x)\n",
    "    for i in range(len(x)):\n",
    "        plt.bar(i, x[i])\n",
    "        plt.xticks( range(5),['other','MY','P','MP','IB'])\n",
    "    mu =np.mean(x)\n",
    "    sigma =np.std(x)\n",
    "    y = norm.pdf(range(15), mu, sigma)#拟合一条最佳正态分布曲线y \n",
    "    print(mu,sigma,y)\n",
    "    # plt.plot(range(15), y, 'r--') #绘制y的曲线\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## mainbranch\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "202568 001.swc\n",
      "210731 043.swc\n",
      "211185 007.swc\n",
      "211185 010.swc\n",
      "211186 001.swc\n",
      "211186 003.swc\n",
      "211188 012.swc\n",
      "211467 016.swc\n",
      "211467 017.swc\n",
      "211467 018.swc\n",
      "211467 019.swc\n",
      "211467 023.swc\n",
      "211467 028.swc\n",
      "211467 032.swc\n",
      "211467 034.swc\n",
      "211468 039.swc\n",
      "211468 062.swc\n",
      "211468 069.swc\n",
      "211470 048.swc\n",
      "211472 004.swc\n",
      "211473 002.swc\n",
      "211988 021.swc\n",
      "211466 016.swc\n",
      "211467 013.swc\n",
      "210728 045.swc\n"
     ]
    }
   ],
   "source": [
    "import sys,copy,os,inspect\n",
    "if hasattr(sys.modules[__name__], '__file__'):\n",
    "    _file_name = __file__\n",
    "else:\n",
    "    _file_name = inspect.getfile(inspect.currentframe())\n",
    "CURRENT_FILE_PATH = os.path.dirname(_file_name)\n",
    "sys.path.append(os.getcwd()+\"/../neuronVis\")\n",
    "\n",
    "import IONData,Scene,SwcLoader\n",
    "import NeuronProcess\n",
    "from sklearn.cluster import KMeans\n",
    "import numpy as np\n",
    "import joblib\n",
    "import pandas as pd\n",
    "\n",
    "def loadClusterCSV(filename,cluster=None):\n",
    "    neurons = pd.read_csv(filename)\n",
    "    neuronsArray= neurons.to_numpy()\n",
    "    neuronScene=[]\n",
    "    for neuron in neuronsArray:\n",
    "        # print(str(neuron[2])[0])\n",
    "        neurondict={}\n",
    "        neurondict['cluster'] = neuron[1]\n",
    "        neurondict['sampleid'] = str(neuron[2])[0:6]\n",
    "        neurondict['name']=str(neuron[2])[6:]+'.swc'\n",
    "        if cluster ==neuron[1] or cluster ==None:\n",
    "            neuronScene.append(neurondict)\n",
    "    return neuronScene\n",
    "\n",
    "def GetFirstOrderEdges(edgesForCluster,parentedge):\n",
    "    for child in parentedge.children:\n",
    "        if child.order!=0:\n",
    "            edge = SwcLoader.Edge()\n",
    "            edgesForCluster.append(edge)\n",
    "            for p in child.data:\n",
    "                edge.addPoint(SwcLoader.Point([p.index,p.type,p.x,p.y,p.z,p.ratio,p.parentIndex]))\n",
    "            child.maxDepth=1\n",
    "            GetFirstOrderEdges(edgesForCluster,child)\n",
    "\n",
    "\n",
    "import Visual as nv\n",
    "import json\n",
    "neuronvis = nv.neuronVis(renderModel=0)\n",
    "\n",
    "neuronvis.render.setBackgroundColor((0.0,0.0,0.,1.0))\n",
    "iondata = IONData.IONData()\n",
    "neurons = loadClusterCSV('../resource/cluster_eachNeuron/cluster_eachNeuron_spcd.csv',cluster=8)\n",
    "neuronvis.render.setLookAt((5000,-15000,0),(5000,0,0),(0,0,-1))\n",
    "mainbranchs=[]\n",
    "edgesForCluster=[]\n",
    "for neuron in neurons[7:-1]:\n",
    "    print(neuron['sampleid'], neuron['name'])\n",
    "    # neuronT = iondata.getNeuronTreeByID('211467','023.swc')\n",
    "    neuronT = iondata.getNeuronTreeByID(neuron['sampleid'], neuron['name'])\n",
    "    neuronT.dendriteHide=True\n",
    "    NeuronProcess.CalculateBranchMaxDepth(neuronT)\n",
    "    mainbranch=[]\n",
    "    if neuronT.rootAxonEdge is None:\n",
    "        continue\n",
    "    for edge in neuronT.edges:\n",
    "\n",
    "        if edge.maxDepth>neuronT.rootAxonEdge.maxDepth*0.5:\n",
    "            edge.order=0\n",
    "            mainbranch.append(edge)\n",
    "    for edge in mainbranch:\n",
    "        for child in edge.children:\n",
    "            if child.order!=0:\n",
    "                NeuronProcess.OrderChildren(neuronT,child,1)\n",
    "        pass\n",
    "    # for child in neuronT.rootAxonEdge.children:\n",
    "    # \t# if child.order!=0:\n",
    "    # \t\tNeuronProcess.OrderChildren(neuronT,child,1)\n",
    "\n",
    "    for edge in mainbranch:\n",
    "        for child in edge.children:\n",
    "            if child.order!=0:\n",
    "                GetFirstOrderEdges(edgesForCluster,child)\n",
    "    mainbranchs.append(mainbranch)\n",
    "    newneuron=SwcLoader.NeuronTree()\n",
    "    if len(mainbranch):\n",
    "        newneuron.edges= mainbranch\n",
    "        newneuron.rootAxonEdge = neuronT.rootAxonEdge\n",
    "        newneuron.root=neuronT.root\n",
    "    neuronvis.addNeuronTree(newneuron,'neuronname',depthIntensity=False)\n",
    "\n",
    "neuronvis.render.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'GeometryAdapter'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32md:\\project\\python\\neuron-vis\\figures\\mainbranch.ipynb Cell 7\u001b[0m in \u001b[0;36m<cell line: 1>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> <a href='vscode-notebook-cell:/d%3A/project/python/neuron-vis/figures/mainbranch.ipynb#W6sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m \u001b[39mimport\u001b[39;00m \u001b[39mGeometryAdapter\u001b[39;00m\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/project/python/neuron-vis/figures/mainbranch.ipynb#W6sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m neuronvis \u001b[39m=\u001b[39m nv\u001b[39m.\u001b[39mneuronVis(renderModel\u001b[39m=\u001b[39m\u001b[39m0\u001b[39m)\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/project/python/neuron-vis/figures/mainbranch.ipynb#W6sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m neuronvis\u001b[39m.\u001b[39mrender\u001b[39m.\u001b[39msetBackgroundColor((\u001b[39m1.0\u001b[39m,\u001b[39m1.0\u001b[39m,\u001b[39m1.\u001b[39m,\u001b[39m1.0\u001b[39m))\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'GeometryAdapter'"
     ]
    }
   ],
   "source": [
    "\n",
    "import GeometryAdapter\n",
    "\n",
    "neuronvis = nv.neuronVis(renderModel=0)\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.,1.0))\n",
    "neuronvis.render.setLookAt((5000,-15000,0),(5000,0,0),(0,0,-1))\n",
    "points=[]\n",
    "branchpoints=[]\n",
    "for mainbranch in mainbranchs[:]:\n",
    "    for edge in mainbranch:\n",
    "        pre = None\n",
    "        for p in edge.data:\n",
    "            if pre is not None :\n",
    "                if (pre.z-5700)*(p.z-5700)<0:\n",
    "                    points.append(p)\n",
    "\n",
    "                    count=0\n",
    "                    parent = edge.parent\n",
    "                    if parent:\n",
    "                        while count<1:\n",
    "                            for child in parent.children:\n",
    "                                if child in mainbranch:\n",
    "                                    count+=1\n",
    "                            if count>1:\n",
    "                                branchpoints.append(parent.data[len(parent.data)-1])\n",
    "                            if parent.parent and count==1:\n",
    "                                parent = parent.parent\n",
    "                                count=0\n",
    "                    \n",
    "                    \n",
    "            pre = p\n",
    "    newneuron.edges= mainbranch\n",
    "    newneuron.rootAxonEdge = neuronT.rootAxonEdge\n",
    "    newneuron.root=neuronT.root\n",
    "    newneuron.somaRadius=100\n",
    "    neuronvis.addNeuronTree(newneuron,name='neuronname',color=[1.0,0.0,0.0],depthIntensity=False)\n",
    "\n",
    "\n",
    "ga = GeometryAdapter.GeometryAdapter()\n",
    "for i in range(len(points)):\n",
    "    ga.geometry.addPoint([points[i].x,points[i].y,points[i].z],[0.0,.0,1.0])\n",
    "    ga.geometry.addIndex(i)\n",
    "ga.geometry.drawModel='points'\n",
    "neuronvis.render.addGeometry(ga.geometry)\n",
    "\n",
    "ga0 = GeometryAdapter.GeometryAdapter()\n",
    "for i in range(len(branchpoints)):\n",
    "    ga0.geometry.addPoint([branchpoints[i].x,branchpoints[i].y,branchpoints[i].z],[0.0,1.0,.0])\n",
    "    ga0.geometry.addIndex(i)\n",
    "ga0.geometry.drawModel='points'\n",
    "neuronvis.render.addGeometry(ga0.geometry)\n",
    "\n",
    "neuronvis.render.setPointSize(10)\n",
    "\n",
    "neuronvis.render.run()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.13 ('neuronVis')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e1da01a652087c3f485a2a8b69806c32fa3b3acf94bcdc0c1ef74e2af5cc794e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
